7  Autocorrelación espacial en ArcGIS y R

Geoestadística i Tècniques d’Observació Global (24/25)

Author

Aitor Ameztegui (UdL)

Published

August 22, 2024

7.1 Introducción

En este tutorial vamos a ver cómo llevar a cabo un análisis de autocorrelación espacial usando tanto R como ArcGIS Pro y QGis.

Esto nos permitirá comparar las diferentes herramientas, analizar en qué difieren y en qué se parecen, y tendremos por tanto un mayor abanico de opciones para realizar este tipo de análisis.

En teoría hemos visto que la autocorrelación espacial no es más que el grado de correlación de una variable consigo misma a través del espacio. Para poder calcularla necesitamos por tanto: 1) observaciones de una variable, y 2) la ubicación espacial de dichas observaciones. Todo ello nos lo proporciona cualquier fichero de información espacial (shapefiles de ESRI, geopackages, u objetos de R de clase sf, como vimos en el tutorial “Accediendo a datos espaciales en R”).

Por ello, lo primero que vamos a hacer es cargar el fichero de datos espaciales que vamos a analizar, que no es otro que el shapefile con los valores de estaciones meteorológicas que usamos en la unidad sobre regresión lineal. También necesitamos cargar el paquete sfdep, que contienen funciones específicas para análisis espaciales en R. Por último, cargaremos el shapefile de provincias para mejorar nuestra visualización:

# Cargamos las librerías que necesitamos
library(sf)
# library(spdep)
library(sfdep) 

# Cargamos el fichero de datos (en este caso, un shapefile)

estaciones <- st_read('data/meteo/meteo_espacial/estaciones_meteo.shp')
Reading layer `estaciones_meteo' from data source 
  `C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\estaciones_meteo.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 90 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N
provincias <- st_read('data/meteo/meteo_espacial/provincias_spain.shp')
Reading layer `provincias_spain' from data source 
  `C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\provincias_spain.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 51 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -14129.47 ymin: 3892590 xmax: 1126923 ymax: 4859517
Projected CRS: ETRS89 / UTM zone 30N

7.2 Definiendo vecinos

Hemos visto en la parte de teoría que uno de los pasos críticos en un análisis de autocorrelación es decidir qué observaciones consideramos como vecinas de un determinado punto. Hay varias manera de definirlo, pero las más importantes son definir los vecinos con criterios de contigüidad (en el caso de polígonos o raster), seleccionar los \(k\) puntos más cercanos, o definir el vecindario en función de la distancia entre los elementos. Métodos más complejos permiten también asignar pesos a los vecinos en función de diversos criterios, entre los que destaca asignarlos de manera inversa a la distancia entre puntos.

7.2.1 Vecinos contiguos

Si tenemos un fichero de polígonos, podemos definir la vecindad como aquellos polígonos contiguos a cada uno, mediante la función st:contiguity() del paquete sfdep. La opción queen = TRUE determina cómo vecinos de un polígono todos aquellos que coincidan en un sólo punto de sus límites. Si queen = FALSE sólo son vecinos aquellos que compartan al menos dos puntos (vecindad de Rook)

veins_cont <- st_contiguity(provincias, queen = FALSE)
veins_cont
Neighbour list object:
Number of regions: 51 
Number of nonzero links: 222 
Percentage nonzero weights: 8.535179 
Average number of links: 4.352941 
4 regions with no links:
7 49 50 51
5 disjoint connected subgraphs

Ahora podemos plotear este objeto para ver las conexiones entre polígonos:

plot(st_geometry(provincias))
plot(veins_cont, coords = st_coordinates(st_centroid(provincias)), add = TRUE, col = "red")
Warning: st_centroid assumes attributes are constant over geometries

NOTA: este método sólo funciona para archivos espaciales de polígonos. Si tratáramos de aplicarlos con puntos obtendríamos un error, ya que los puntos no pueden tener puntos de contacto entre ellos.

NOTA (2): aunque la opción “queen” permite definir el criterio de selección de vecinos, en el caso de polígonos complejos como el del fichero provincias lo más normal es que no haya diferencias entre queen = TRUE o queen = FALSE

7.2.2 \(k\) vecinos más próximos

Un método alternativo consiste en elegir como vecinos los \(k\) puntos más cercanos, lo que se adapta a toda la zona de estudio, teniendo en cuenta las diferencias en las densidades de las entidades de área. En este caso tan sólo debemos definir el valor de \(k\), y la función buscará automáticamente los k puntos más próximos de cada punto. Usaremos para ello la función st_knn():

veins_k <- st_knn(estaciones, k = 8)
veins_k
Neighbour list object:
Number of regions: 90 
Number of nonzero links: 720 
Percentage nonzero weights: 8.888889 
Average number of links: 8 
Non-symmetric neighbours list

Ahora podemos plotearlo como antes:

plot(veins_k, coords = st_coordinates(estaciones), col = "red")

En el caso de esta metodología, puede utilizarse también con polígonos. Sin embargo requiere que trabajemos con los centroides de cada polígono, que serán los que nos permitan determinar la distancia entre polígonos:

veins_prov_k <- st_knn(st_centroid(provincias), k = 4)
Warning: st_centroid assumes attributes are constant over geometries
plot(st_geometry(provincias))
plot(veins_prov_k, coords = st_coordinates(st_centroid(provincias)),
     col = "red", add = TRUE)
Warning: st_centroid assumes attributes are constant over geometries

7.2.3 Vecinos basados en intervalos de distancia

Para definir los vecinos en base a una distancia podemos usar la función st_dist_band() en la que debemos definir la distancia mínima (lower) y máxima (upper) para seleccionar un vecino. La función determinará que todos los puntos que estén entre estas dos distancias serán vecinos:

veins_dist <- st_dist_band(estaciones, lower = 0, upper = 30000)
plot(veins_dist, coords = st_coordinates(estaciones), col = "red")

Igual que antes, si queremos usar este método con polígonos, debemos hacerlo con los centroides de los polígonos.

7.3 Asignando pesos a los vecinos

Otra opción muy interesante, y que suele aplicarse en muchos casos, es no sólo definir cuáles consideraremos como vecinos de cada observación, sino darles un peso determinado según la distancia a la que estén del punto analizado. Este método se conoce como weighting.

Aunque existen varias opciones para dar pesos, veremos las dos más comunes: asignar el mismo peso a todos los vecinos de un punto (función st_weights()) o asignar mayor peso a los vecinos que estén más cerca, lo que se conoce como inverse distance weighting (idw) (función st_inverse_distance()).

Para evitar que la función calcule las distancias entre todos los puntos dos a dos - lo que podría colapsar la memoria del ordenador - las dos funciones se aplican sobre un objeto de vecinos ya creado, en el que le digamos cuáles son los vecinos de cada punto. En este caso, usaremos el que acabamos de crear (veins_dist):

veins_w_dist <- st_weights(veins_dist)
veins_w_dist
[[1]]
[1] 0.5 0.5

[[2]]
[1] 1

[[3]]
[1] 0.2 0.2 0.2 0.2 0.2

[[4]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[5]]
[1] 0.25 0.25 0.25 0.25

[[6]]
[1] 0.3333333 0.3333333 0.3333333

[[7]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[8]]
[1] 0.25 0.25 0.25 0.25

[[9]]
[1] 0.25 0.25 0.25 0.25

[[10]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[11]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[12]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[13]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[14]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[15]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[16]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[17]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[18]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[19]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[20]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[21]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[22]]
 [1] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333
 [7] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333

[[23]]
 [1] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333
 [7] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333

[[24]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[25]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[26]]
 [1] 0.09090909 0.09090909 0.09090909 0.09090909 0.09090909 0.09090909
 [7] 0.09090909 0.09090909 0.09090909 0.09090909 0.09090909

[[27]]
 [1] 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308
 [7] 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308 0.07692308
[13] 0.07692308

[[28]]
[1] 0.2 0.2 0.2 0.2 0.2

[[29]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[30]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[31]]
[1] 0.2 0.2 0.2 0.2 0.2

[[32]]
[1] 0.2 0.2 0.2 0.2 0.2

[[33]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[34]]
 [1] 0.09090909 0.09090909 0.09090909 0.09090909 0.09090909 0.09090909
 [7] 0.09090909 0.09090909 0.09090909 0.09090909 0.09090909

[[35]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[36]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[37]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[38]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[39]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[40]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[41]]
[1] 0.2 0.2 0.2 0.2 0.2

[[42]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[43]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[44]]
[1] 0.25 0.25 0.25 0.25

[[45]]
[1] 0.2 0.2 0.2 0.2 0.2

[[46]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[47]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[48]]
[1] 0.25 0.25 0.25 0.25

[[49]]
[1] 0.25 0.25 0.25 0.25

[[50]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[51]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[52]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[53]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[54]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[55]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[56]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[57]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[58]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[59]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[60]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[61]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[62]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[63]]
[1] 1

[[64]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[65]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[66]]
 [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

[[67]]
[1] 0.25 0.25 0.25 0.25

[[68]]
[1] 0.25 0.25 0.25 0.25

[[69]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[70]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[71]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[72]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[73]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[74]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[75]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[76]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[77]]
[1] 0.2 0.2 0.2 0.2 0.2

[[78]]
[1] 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125

[[79]]
[1] 0.25 0.25 0.25 0.25

[[80]]
[1] 0.3333333 0.3333333 0.3333333

[[81]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

[[82]]
[1] 0.3333333 0.3333333 0.3333333

[[83]]
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

[[84]]
[1] 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111 0.1111111
[8] 0.1111111 0.1111111

[[85]]
[1] 0.25 0.25 0.25 0.25

[[86]]
[1] 0.3333333 0.3333333 0.3333333

[[87]]
[1] 0.3333333 0.3333333 0.3333333

[[88]]
[1] 0.25 0.25 0.25 0.25

[[89]]
[1] 0.3333333 0.3333333 0.3333333

[[90]]
[1] 0.5 0.5

attr(,"mode")
[1] "binary"
attr(,"W")
[1] TRUE
attr(,"comp")
attr(,"comp")$d
 [1]  2  1  5  6  4  3  7  4  4  7  7  8  8  6  8 10  6  9  7  9 10 12 12 10  7
[26] 11 13  5  9  9  5  5  9 11 10  7  6  7  7  9  5  7  9  4  5 10  7  4  4  8
[51]  8  9  9  9  9  7  9  8 10  7  7  6  1  7  6 10  4  4  8  8  6  7  6  7  8
[76]  6  5  8  4  3  7  3  6  9  4  3  3  4  3  2
veins_idw <- st_inverse_distance(veins_dist, estaciones)
veins_idw
[[1]]
[1] 0.006162935 0.004629866

[[2]]
[1] 0.004474438

[[3]]
[1] 0.003614838 0.008617471 0.005888738 0.004308553 0.004246341

[[4]]
[1] 0.004110196 0.005004989 0.005892665 0.005059494 0.003855733 0.003525225

[[5]]
[1] 0.006162935 0.004110196 0.018573094 0.004219539

[[6]]
[1] 0.003614838 0.003762743 0.003742014

[[7]]
[1] 0.010436978 0.003739723 0.007476455 0.005599893 0.003835977 0.004550622
[7] 0.004013768

[[8]]
[1] 0.004629866 0.005004989 0.018573094 0.005150181

[[9]]
[1] 0.004474438 0.003358967 0.003759661 0.004728103

[[10]]
[1] 0.008617471 0.003762743 0.012830421 0.006860108 0.008371453 0.004122217
[7] 0.003432164

[[11]]
[1] 0.005892665 0.035738305 0.006207863 0.003753747 0.007777664 0.003519870
[7] 0.003568408

[[12]]
[1] 0.010436978 0.003358967 0.011589005 0.012081789 0.005989138 0.005008006
[7] 0.006173371 0.003603735

[[13]]
[1] 0.005059494 0.035738305 0.007278701 0.003583814 0.009221188 0.003575821
[7] 0.003899278 0.003683868

[[14]]
[1] 0.005888738 0.012830421 0.014673056 0.009628569 0.004394696 0.004420994

[[15]]
[1] 0.003739723 0.006207863 0.007278701 0.003904562 0.006026709 0.005189466
[7] 0.005225900 0.005394335

[[16]]
 [1] 0.007476455 0.011589005 0.003904562 0.008790135 0.004959863 0.008764587
 [7] 0.007545888 0.004368521 0.003813779 0.004515097

[[17]]
[1] 0.004308553 0.006860108 0.014673056 0.007949523 0.004408820 0.005841640

[[18]]
[1] 0.005599893 0.003759661 0.012081789 0.008790135 0.011061750 0.004623875
[7] 0.010092823 0.004479795 0.003718229

[[19]]
[1] 0.004246341 0.008371453 0.009628569 0.007949523 0.003721542 0.007871388
[7] 0.005058706

[[20]]
[1] 0.003835977 0.004728103 0.005989138 0.004959863 0.011061750 0.003403491
[7] 0.007795531 0.004218550 0.004983057

[[21]]
 [1] 0.003855733 0.004219539 0.005150181 0.003753747 0.003583814 0.003795294
 [7] 0.003713069 0.004750135 0.003414890 0.004118014

[[22]]
 [1] 0.003525225 0.007777664 0.009221188 0.006026709 0.003795294 0.003402316
 [7] 0.004575316 0.005373133 0.005561090 0.004853792 0.004557856 0.003879382

[[23]]
 [1] 0.004550622 0.005008006 0.005189466 0.008764587 0.004623875 0.003403491
 [7] 0.003402316 0.005418191 0.008705479 0.006739940 0.004848548 0.003681319

[[24]]
 [1] 0.004013768 0.006173371 0.007545888 0.010092823 0.007795531 0.005418191
 [7] 0.003636457 0.008045509 0.004538472 0.003611401

[[25]]
[1] 0.003742014 0.003721542 0.005031876 0.004780327 0.003337743 0.004294985
[7] 0.003475826

[[26]]
 [1] 0.003575821 0.005225900 0.004368521 0.004575316 0.008705479 0.003636457
 [7] 0.026493047 0.004075294 0.003516125 0.003962043 0.004654760

[[27]]
 [1] 0.003519870 0.003899278 0.005394335 0.003813779 0.005373133 0.006739940
 [7] 0.026493047 0.003516832 0.003571777 0.004019564 0.004558003 0.003638788
[13] 0.004400688

[[28]]
[1] 0.023162289 0.003507235 0.004820755 0.003749024 0.004364900

[[29]]
[1] 0.004122217 0.004394696 0.004408820 0.007871388 0.005031876 0.005210361
[7] 0.004499440 0.004485682 0.003400160

[[30]]
[1] 0.003432164 0.004420994 0.005841640 0.005058706 0.005210361 0.003436339
[7] 0.005719380 0.003931030 0.003448073

[[31]]
[1] 0.003713069 0.007232551 0.006956508 0.005500667 0.004358544

[[32]]
[1] 0.023162289 0.003505615 0.005370645 0.004124959 0.005333420

[[33]]
[1] 0.003568408 0.003683868 0.004750135 0.005561090 0.003516832 0.012147185
[7] 0.003587408 0.007556955 0.007787294

[[34]]
 [1] 0.003603735 0.004515097 0.004479795 0.004218550 0.004848548 0.008045509
 [7] 0.004075294 0.003571777 0.004457251 0.005070455 0.004870103

[[35]]
 [1] 0.003718229 0.004983057 0.004538472 0.003507235 0.003505615 0.004457251
 [7] 0.006659033 0.006156431 0.006765926 0.003536958

[[36]]
[1] 0.003414890 0.004853792 0.003516125 0.004019564 0.012147185 0.019100230
[7] 0.018443203

[[37]]
[1] 0.004118014 0.007232551 0.003587408 0.004689643 0.011410470 0.004002355

[[38]]
[1] 0.004557856 0.003962043 0.004558003 0.007556955 0.019100230 0.018040963
[7] 0.003398533

[[39]]
[1] 0.003436339 0.006956508 0.004689643 0.005252977 0.009941154 0.004210954
[7] 0.005292588

[[40]]
[1] 0.004820755 0.005370645 0.006659033 0.016418318 0.004811032 0.006453718
[7] 0.004859959 0.004620294 0.003536070

[[41]]
[1] 0.003879382 0.003638788 0.007787294 0.018443203 0.018040963

[[42]]
[1] 0.004780327 0.004499440 0.003967366 0.003735427 0.010401674 0.012143919
[7] 0.003705614

[[43]]
[1] 0.004485682 0.005719380 0.003967366 0.003765583 0.006117359 0.004333852
[7] 0.004748846 0.004193798 0.003862178

[[44]]
[1] 0.005500667 0.011410470 0.005252977 0.005241227

[[45]]
[1] 0.003681319 0.004654760 0.004400688 0.005070455 0.003398533

[[46]]
 [1] 0.003749024 0.004124959 0.006156431 0.016418318 0.005868575 0.006075711
 [7] 0.006788932 0.005458624 0.004345378 0.003370678

[[47]]
[1] 0.003611401 0.004870103 0.006765926 0.004811032 0.005868575 0.004925559
[7] 0.003394432

[[48]]
[1] 0.003337743 0.003735427 0.005663043 0.003707812

[[49]]
[1] 0.004294985 0.010401674 0.005663043 0.010005842

[[50]]
[1] 0.004364900 0.005333420 0.006453718 0.006075711 0.004509175 0.007420818
[7] 0.004023550 0.004513244

[[51]]
[1] 0.004358544 0.004002355 0.009941154 0.005241227 0.005155981 0.007888297
[7] 0.003451690 0.003391276

[[52]]
[1] 0.003475826 0.003400160 0.012143919 0.003765583 0.003707812 0.010005842
[7] 0.004611517 0.003773790 0.003696297

[[53]]
[1] 0.003931030 0.004210954 0.006117359 0.005155981 0.013964334 0.004101068
[7] 0.007014490 0.005723988 0.005219137

[[54]]
[1] 0.003448073 0.005292588 0.004333852 0.007888297 0.013964334 0.005808662
[7] 0.004894253 0.004986876 0.003720591

[[55]]
[1] 0.003536958 0.004859959 0.006788932 0.004925559 0.004509175 0.007608761
[7] 0.010541089 0.004673287 0.004478467

[[56]]
[1] 0.004620294 0.005458624 0.007420818 0.007608761 0.008415016 0.008783083
[7] 0.003614692

[[57]]
[1] 0.003705614 0.004748846 0.004611517 0.004101068 0.005811761 0.006178186
[7] 0.004559716 0.004588461 0.004121207

[[58]]
[1] 0.004193798 0.003451690 0.007014490 0.005808662 0.005811761 0.030537597
[7] 0.015206835 0.004123002

[[59]]
 [1] 0.003536070 0.004345378 0.003394432 0.004023550 0.010541089 0.008415016
 [7] 0.006892460 0.004938095 0.004549184 0.004713823

[[60]]
[1] 0.003862178 0.005723988 0.004894253 0.006178186 0.030537597 0.017182885
[7] 0.004091558

[[61]]
[1] 0.003370678 0.004513244 0.004673287 0.008783083 0.006892460 0.005466388
[7] 0.004309843

[[62]]
[1] 0.005219137 0.004986876 0.004559716 0.015206835 0.017182885 0.005361710

[[63]]
[1] 0.004029734

[[64]]
[1] 0.003500060 0.003645244 0.005738888 0.003351490 0.005768619 0.003500076
[7] 0.003407870

[[65]]
[1] 0.004478467 0.004938095 0.004029734 0.003357246 0.004378309 0.003531956

[[66]]
 [1] 0.003391276 0.003720591 0.004123002 0.004091558 0.005361710 0.005095665
 [7] 0.004811194 0.007141728 0.004621374 0.003354525

[[67]]
[1] 0.003773790 0.004588461 0.035889025 0.003398814

[[68]]
[1] 0.003696297 0.004121207 0.035889025 0.003489890

[[69]]
[1] 0.003500060 0.005095665 0.079968735 0.010158248 0.027584200 0.008114575
[7] 0.007090058 0.003383981

[[70]]
[1] 0.003645244 0.004811194 0.079968735 0.009066521 0.027098256 0.008602432
[7] 0.007182009 0.003501714

[[71]]
[1] 0.005738888 0.011110630 0.005299489 0.005785850 0.004056352 0.003574184

[[72]]
[1] 0.003614692 0.004549184 0.005466388 0.003357246 0.011662586 0.007968841
[7] 0.004110888

[[73]]
[1] 0.007141728 0.010158248 0.009066521 0.010301497 0.005602327 0.006286845

[[74]]
[1] 0.004713823 0.004309843 0.004378309 0.011662586 0.010048597 0.003970192
[7] 0.003596182

[[75]]
[1] 0.003351490 0.004621374 0.027584200 0.027098256 0.010301497 0.010704627
[7] 0.009532459 0.003679035

[[76]]
[1] 0.005768619 0.011110630 0.004069532 0.003852847 0.005600518 0.006262039

[[77]]
[1] 0.005299489 0.004069532 0.007636479 0.007900458 0.009344532

[[78]]
[1] 0.003500076 0.008114575 0.008602432 0.005602327 0.010704627 0.003852847
[7] 0.016374081 0.005518690

[[79]]
[1] 0.007968841 0.010048597 0.004385258 0.005127903

[[80]]
[1] 0.010997655 0.003568856 0.004974987

[[81]]
[1] 0.003354525 0.007090058 0.007182009 0.006286845 0.009532459 0.016374081
[7] 0.004673300

[[82]]
[1] 0.007636479 0.004542729 0.011909997

[[83]]
[1] 0.005785850 0.005600518 0.007900458 0.004542729 0.004200314 0.006786618

[[84]]
[1] 0.003407870 0.003383981 0.003501714 0.004056352 0.003679035 0.006262039
[7] 0.005518690 0.004673300 0.004200314

[[85]]
[1] 0.003574184 0.009344532 0.011909997 0.006786618

[[86]]
[1] 0.003531956 0.003970192 0.004385258

[[87]]
[1] 0.010997655 0.004585941 0.004477688

[[88]]
[1] 0.003398814 0.003489890 0.003568856 0.004585941

[[89]]
[1] 0.004110888 0.003596182 0.005127903

[[90]]
[1] 0.004974987 0.004477688

Al imprimir el resultado vemos que, mientras en veins_dist_w() todos los vecinos de una observación tienen el mismo peso, en veins_idw los pesos varían entre los vecinos.

7.4 Autocorrelació Global: I de Moran

La I global de Moran determina el grado de autocorrelación de una muestra en su conjunto. Valores positivos indican autocorrelación positiva (los valores más próximos entre sí muestran valores más similares), valores negativos muestran lo contrario, y valores cercanos a 0, ausencia de autocorrelación espacial.

7.4.1 I global de Moran con R

En R podemos determinar la I de Moran con la función global_moran_test() que necesita como argumentos: (1) la variable que queremos analizar (en este caso será la temperatura media del mes, Tmed_MES); (2) un objeto neighbor que identifique los vecinos de cada observación; y (3) un objeto de pesos. Por lo tanto los pasos a seguir son:

  1. Identificar los vecinos de cada observación
  2. Asignar los pesos mediante la opción st_weights() o st_inverse_distance()
  3. Calcular la I de Moran

Aunque ya tenemos los objetos resultantes de los pasos 1 y 2, vamos a repasar todo el proceso para tener todos los pasos juntos:

7.4.1.1 I global de Moran en base a los k vecinos más próximos

Como hemos visto antes, en el caso de usar los k vecinos más próximos, usaremos las funciones st_knn() y st_weights():

# 1. Identificar los vecinos de cada observación en un objeto neighbour (nb)
veins_k <- st_knn(estaciones, k = 8)

# 2. Asignar los pesos de cada observación (en este caso, pesos iguales)
veins_w_k <- st_weights(veins_k)

# 3. Calcular la I de Moran para la variable Tmed_MES
globalMoran_k <- global_moran_test(x = estaciones$Tmed_MES, nb = veins_k, wt = veins_w_k)
globalMoran_k

    Moran I test under randomisation

data:  x  
weights: listw    

Moran I statistic standard deviate = 12.934, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.607761749      -0.011235955       0.002290385 

Vemos que el valor de la I de Moran es 0.6077, lo que indica autocorrelación espacial positiva. El valor esperado si los datos se distribuyeran aleatoriamente sería de -0.011. Por último, también sabemos que la probabilidad de que los valores de temperatura observados estén distribuidos al azar es de 2.2e-16, con lo que claramente rechazaremos la hipótesis nula (es decir, rechazamos que no haya autocorrelación). En definitiva, que sí que hay autocorrelación espacial.

7.4.1.2 I global de Moran en base a la distancia entre vecinos

En este caso reconoceremos como vecinos aquellos dentro de un intervalo de distancias, y asignaremos el mismo peso a todos:

# 1. Identificar los vecinos de cada observación en un objeto neighbour (nb)
veins_dist <- st_dist_band(estaciones, lower = 0, upper = 30000)

# 2. Asignar los pesos de cada observación (en este caso, pesos iguales)
veins_w_dist <- st_weights(veins_dist)

# 3. Calcular la I de Moran
globalMoran_dist <- global_moran_test(x = estaciones$Tmed_MES, nb = veins_dist, wt = veins_w_dist)
globalMoran_dist

    Moran I test under randomisation

data:  x  
weights: listw    

Moran I statistic standard deviate = 11.134, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.664310960      -0.011235955       0.003681663 

Vemos que el valor cambia respecto al calculado antes (0.664),ya que la definición de cuáles son los vecinos ha cambiado mucho. Sin embargo, la conclusión sigue siendo la misma (hay bastante autocorrelación)

7.4.1.3 I Global de Moran en base a IDW

Ahora identificaremos los vecinos en base a la franja de distancias, pero les asignaremos un peso inversamente proporcional a la distancia con el punto analizado:

# 1. Identificar los vecinos de cada observación en un objeto neighbour (nb)
veins_dist <- st_dist_band(estaciones, lower = 0, upper = 30000)

# 2. Asignar los pesos de cada observación (en este caso, pesos iguales)
veins_idw <- st_inverse_distance(veins_dist, estaciones)

# 3. Calcular la I de Moran
globalMoran_idw <- global_moran_test(x = estaciones$Tmed_MES, nb = veins_dist, wt = veins_idw)
globalMoran_idw

    Moran I test under randomisation

data:  x  
weights: listw    

Moran I statistic standard deviate = 8.7647, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.663803406      -0.011235955       0.005931776 

Aunque el resultado es parecido al de globalMoran_dist, hay pequeñas diferencias debidas a los pesos asignados. En función del dataset concreto las diferencias pueden ser más o menos grandes.

7.4.2 I Global de Mo#ran con ArcGis Pro

Determinar la I global de Moran con ArcGIS es también muy sencillo. En primer lugar, debemos cargar la capa estaciones_meteo.shp, y tenemos que buscar la función Spatial Autocorrelation (Global Moran’s I), que podemos encontrar dentro de las herramientas de Geoprocessing

Una vez seleccionada la función, tenemos que indicar la capa que queremos analizar, la variable de interés, y cómo determinaremos el peso de los vecinos.

Los resultados se muestran en la pestaña Results de ArcGIS, que si no tenemos visible, podemos activar yendo a Geoprocessing/Results. Si además hemos seleccionado que queremos que genere un informe, este estará disponible como html, produciendo una salida visualmente muy informativa.

> ** Ejercicio** Prueba a modificar la manera de definir los pesos de cada observación, verás que los resultados de ArcGIS son los mismos que obtenemos en R.

7.4.3 I Global de Moran con QGis y Geoda

7.4.3.1 Presentando GeoDa

Por supuesto, QGis tiene funcionalidades parecidas a las de ArcGis, mediante un plugin llamado Spatial Analysis Toolbox. Sin embargo, a día de hoy (julio de 2024), este plugin no funciona, así que veremos una alternativa muy interesante, el software GeoDa.

GeoDa es un software de escritorio gratuito y de código abierto, desarrollado por el equipo del Dr. Luc Anselin, profesor de la Universidad de Chicago y uno de los grandes pioneros en el análisis espacial y la econometría. Podemos descargar gratuitamente GeoDa aquí.

Una vez descargado, si lo abrimos nos pedirá que carguemos los datos: podemos cargar capas y objetos espaciales en multitud de formatos. En este caso, vamos a cargar el fichero estaciones_meteo.shp.

Una vez cargados los datos, el aspecto es el de un SIG convencional, con la tabla de contenidos a la izquierda, la visualización a la derecha, y arriba veremos el menú con diversas opciones.

7.4.3.2 Calculando vecinos y pesos en GeoDa

Para poder calcular la I de Moran, GeoDa nos exige que calculemos los pesos. Para ello debemos acceder al menú Weights Manager caracterizado por una gran W morada. El menú del Weights manager nos resultará familiar, ya que nos da la opción de definir los vecinos por contiguidad (según el criterio de Rooke o de Queen) o por distancia (y en este caso podremos elegir por banda de distancia, k vecinos más próximos o mediante un kernel, que no lo hemos visto). En caso de elegir los vecinos por distancia, nos preguntará si queremos usar el inverso de la distancia como peso (es decir, un idw):

Una vez decidido el criterio que queramos, guardamos el fichero de pesos en disco, con el nombre que queramos y la extensión .gwd. Una de las funcionalidades más interesantes de GeoDa es que nos permite visualizar, una vez seleccionado el fichero de pesos que queramos, los vecinos de un punto sobre el que pongamos el cursos, además de permitirnos ver gráficos de vecindad interactivos:

Podemos probar a crear un fichero con knn (8 vecinos), otro con los vecinos entre 0 y 30.000 metros, y otro con los mismos vecinos pero pesos inversos a la distancia, como hicimos en R.

7.4.3.3 Calculando la I de Moran con GeoDa

Una vez definidos los vecinos, podemos calcular fácilmente la I global de Moran, seleccionando el comando Univariate Moran's I del menú Space. Para ello sólo debemos indicar la variable de interés (Tmed_MES) y el fichero de vecinos deseado:

7.5 Autocorrelación Local (I de Anselin)

Para poder calcular la I local de Moran (también llamada I de Anselin), necesitaremos, igual que antes, conocer la variable a utilizar, la identificación de los vecinos, y los pesos. Vamos a trabajar, por simplicidad, con uno de los casos (idw), pero funcionaría igual con el resto:

localMoran <- local_moran(x = estaciones$Tmed_MES, nb = veins_dist, wt = veins_idw)

Vemos que en este caso el objeto generado contiene el valor calculado de I (ii), el valor esperado (eii), la varianza (var_ii), el z-score (z_ii) y el p-valor (p_ii) de cada observación de la muestra, entre otros campos

head(localMoran)
          ii           eii       var_ii     z_ii         p_ii p_ii_sim
1 0.02504040 -1.240126e-03 2.145584e-04 1.794161 0.0727875371    0.096
2 0.01013481  4.419634e-05 3.629683e-05 1.674879 0.0939579963    0.100
3 0.07910570 -1.612636e-03 5.279889e-04 3.512852 0.0004433249    0.004
4 0.03636902  5.256534e-04 2.594232e-04 2.225381 0.0260556745    0.032
5 0.05093252  2.016048e-05 5.892954e-04 2.097281 0.0359686814    0.032
6 0.03242147  1.667611e-04 1.122603e-04 3.044248 0.0023326282    0.004
  p_folded_sim    skewness    kurtosis     mean   median     pysal
1        0.048  0.22413422 -0.45466255  Low-Low  Low-Low   Low-Low
2        0.054 -0.23127189 -0.91338204 High-Low High-Low High-High
3        0.002  0.25155287 -0.08314579  Low-Low  Low-Low   Low-Low
4        0.016  0.10151188 -0.06217115  Low-Low  Low-Low   Low-Low
5        0.016  0.20330568 -0.73921216  Low-Low  Low-Low   Low-Low
6        0.002 -0.00244472 -0.61948853  Low-Low  Low-Low   Low-Low

7.5.1 LISA (Local indicators of spatial association) - Agrupaciones espaciales

Un análisis muy interesante cuando trabajamos con autocorrelación espacial local es el de identificar aquellas observaciones en las que se cumplen una serie de condiciones. Este análisis fue desarrollado en 1995 por Luc Anselin, y recibe el nombre de LISA (local indicators of spatial association). Se compone de varios elementos:

7.5.1.1 Scatterplot de Anselin

Si queremos producir un scatterplot con los valores positivamente y negativamente autocorrelacionados, podemos usar la función moran.plot(), del paquete spdep. Por desgracia, aún no existe una versión de esta función para sfdep, por lo que nos vemos obligados a definir los pesos de los vecinos usando la función nb2listw() (en caso de querer aplicar pesos iguales), o la función nb2listwdist() si queremos aplicar un idw.

library(spdep)
Cargando paquete requerido: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
moran.plot(x = estaciones$Tmed_MES, nb2listw(veins_dist))

moran.plot(x = estaciones$Tmed_MES, nb2listwdist(veins_dist, estaciones))

Vemos que los valores de las esquinas inferior izquierda y superior derecha son los que tienen autocorrelación positiva, y los de las esquinas contrarias, autocorrelación negativa.

7.5.1.2 Mapa de valores de I

También podemos usar la función tm_shape() y tm_dots() del paquete tmap para visualizar espacialmente los valores locales de I de Moran. Primero debemos unir los cálculos de I local al objeto espacial que tenemos, usando cbind().

lisa <- cbind(estaciones, localMoran)

Después podemos crear la representación espacial de la siguiente manera:

library(tmap)
Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
tm_shape(lisa) +
    tm_dots(col = "ii", size = 0.75) +
tm_shape(provincias) +
    tm_borders()
Variable(s) "ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

7.5.1.3 Mapa de agrupaciones cluster

Consiste en identificar las observaciones que cumplen estos criterios:

  • Alto valor de la variable y vecinos con alto valor (High-high)

  • Bajo valor de la variable y vecinos con bajo valor (Low-Low)

  • Alto valor de la variable y vecinos con bajo valor (High-Low)

  • Bajo valor de la variable y vecinos con alto valor (Low-High)

Los dos primeros se suelen colorear de color rojo y azul intenso, respectivamente, y corresponden a las observaciones que contribuyen significativamente a la autocorrelación espacial global positiva. Los otros dos se suelen representar con colores pálidos, y son las observaciones que contribuyen a una autocorrelación negativa. Las observaciones con valores no significativos no suelen colorearse.

La función localMoran ya nos ofrece una clasificación en estos cuatro grupos, a la que podemos acceder de la siguiente manera:

localMoran$mean
 [1] Low-Low   High-Low  Low-Low   Low-Low   Low-Low   Low-Low   Low-Low  
 [8] Low-Low   High-Low  Low-Low   Low-High  Low-High  Low-High  Low-Low  
[15] High-Low  High-High Low-Low   Low-High  Low-Low   High-High Low-Low  
[22] Low-High  High-High High-High Low-Low   High-High High-High High-Low 
[29] Low-Low   Low-Low   Low-Low   High-High High-High High-High High-High
[36] High-High Low-Low   High-High Low-Low   High-High High-High Low-Low  
[43] Low-Low   Low-Low   High-Low  High-High High-Low  Low-Low   Low-Low  
[50] High-High Low-Low   Low-Low   Low-High  Low-High  High-High Low-High 
[57] Low-Low   Low-High  High-High High-High High-Low  High-High High-Low 
[64] High-Low  High-Low  Low-High  Low-Low   Low-Low   High-High High-High
[71] High-Low  High-Low  High-High High-High High-High High-Low  High-Low 
[78] High-High High-Low  High-Low  High-High High-Low  High-Low  High-Low 
[85] High-Low  High-Low  High-Low  Low-Low   High-Low  High-Low 
Levels: Low-Low High-Low Low-High High-High

Como antes hemos unido los resultados de localMoran a nuestra tabla estaciones también tendremos esta información en el objeto moran.map

Si definimos un nivel de significación, podemos clasificar las observaciones según el grupo al que pertenezcan y su significación. Indiquemos que los valores no significativos tengan un valor de “NS”

# llindar de significacion
signif <- 0.1 

lisa$mean <- ifelse(lisa$p_ii > signif,  "NS", as.character(lisa$mean))
lisa
Simple feature collection with 90 features and 25 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
   OBJECTID FID_weathe INDICATIVO AÃ_O MES                    NOMBRE ALTITUD
1         1         12      9382E 2012   6         PANCRUDO (D.G.A.)    1285
2         2         85      9987N 2012   6   DELTEBRE (PARC NATURAL)       4
3         3         11      9377E 2012   6               EL PEDREGAL     985
4         4         50      9567U 2012   6           EJULVE (D.G.A.)    1095
5         5         39      9530V 2012   6         UTRILLAS (D.G.A.)     960
6         6          0       3013 2012   6          MOLINA DE ARAGON    1063
7         7         70      9939A 2012   6        FUENTESPALDA (DGA)     720
8         8         40      9531X 2012   6    MONTALBAN 'AUTOMATICA'     885
9         9         84      9981A 2012   6 TORTOSA (OBSER. DEL EBRO)      48
10       10         89       9999 2012   6                      ODON    1110
   T_MAX_abs T_MIN_abs TMed_MAX TMed_MIN Tmed_MES   Provincia           ii
1       29.5       3.0     22.0      9.6    15.80      Teruel 0.0250404031
2       32.0      10.0     27.5     17.4    22.45   Tarragona 0.0101348145
3       31.0       3.5     21.5     10.2    15.85 Guadalajara 0.0791057035
4       30.0       5.0     21.9     11.5    16.70      Teruel 0.0363690201
5       32.5       3.5     23.7     10.8    17.25      Teruel 0.0509325181
6       32.8       3.2     23.8      8.4    16.10 Guadalajara 0.0324214727
7       32.0       5.0     24.4     13.2    18.80      Teruel 0.0006529013
8       31.2       5.4     23.0     11.0    17.00      Teruel 0.0515590522
9       34.4      13.0     29.3     17.1    23.20   Tarragona 0.0048631840
10      32.0       4.0     22.4     10.2    16.30      Teruel 0.1221026551
             eii       var_ii      z_ii         p_ii p_ii_sim p_folded_sim
1  -1.240126e-03 2.145584e-04 1.7941606 0.0727875371    0.096        0.048
2   4.419634e-05 3.629683e-05 1.6748789 0.0939579963    0.100        0.054
3  -1.612636e-03 5.279889e-04 3.5128517 0.0004433249    0.004        0.002
4   5.256534e-04 2.594232e-04 2.2253810 0.0260556745    0.032        0.016
5   2.016048e-05 5.892954e-04 2.0972813 0.0359686814    0.032        0.016
6   1.667611e-04 1.122603e-04 3.0442479 0.0023326282    0.004        0.002
7  -1.191728e-04 4.568722e-05 0.1142251 0.9090593525    0.916        0.458
8  -1.541596e-04 7.247982e-04 1.9208478 0.0547508938    0.068        0.034
9  -9.059528e-04 1.992059e-04 0.4087519 0.6827217648    0.716        0.358
10 -1.976345e-03 1.113684e-03 3.7180670 0.0002007530    0.004        0.002
      skewness    kurtosis     mean   median     pysal               geometry
1   0.22413422 -0.45466255  Low-Low  Low-Low   Low-Low POINT (666121 4514365)
2  -0.23127189 -0.91338204 High-Low High-Low High-High POINT (814489 4514857)
3   0.25155287 -0.08314579  Low-Low  Low-Low   Low-Low POINT (620767 4515276)
4   0.10151188 -0.06217115  Low-Low  Low-Low   Low-Low POINT (705449 4516865)
5   0.20330568 -0.73921216  Low-Low  Low-Low   Low-Low POINT (681326 4520030)
6  -0.00244472 -0.61948853  Low-Low  Low-Low   Low-Low POINT (593975 4522166)
7   0.17087867 -0.31646450       NS  Low-Low   Low-Low POINT (758868 4522277)
8   0.27739685 -0.58833401  Low-Low  Low-Low   Low-Low POINT (686217 4522281)
9  -0.20533533 -0.18335013       NS High-Low High-High POINT (794466 4524785)
10  0.17227708 -0.13938989  Low-Low  Low-Low   Low-Low POINT (620133 4526863)

Por último, ploteamos la figura con tmap, tm_shape() y tm_dots():

# plot dels resultats

tm_shape(lisa) +
    tm_dots(col = "mean",
          style = "cat",
          palette= c("NS"="grey",
                     "Low-Low" = "blue",        # Low-Low
                     "Low-High" = "lightblue",  # Low-High 
                     "High-Low"= "sienna1",     # High-low
                     "High-High"="red"),        # High -High
              size = 0.5) +
tm_shape(provincias)+
    tm_borders()

7.5.1.4 Mapa de significación

Otro elemento interesante es producir un mapa de significación. Para ello establecemos las categorías de significación que queramos, y ploteamos su distribución espacial. Por ejemplo, vamos a crear 4 categorías:

  • No significativo (P > 0.05)

  • 0.01 < P < 0.05

  • 0.001 < P < 0.01

  • P < 0.001

lisa$sign <- ifelse(lisa$p_ii < 0.001, "p<0.001",
                         ifelse(lisa$p_ii < 0.01, "p<0.01",
                                ifelse(lisa$p_ii < 0.05, "p<0.05",
                                "NS")))

Y creamos el mapa de manera muy similar:

tm_shape(lisa) +
    tm_dots(col = "sign",
          style = "cat",
          palette= c("NS"="grey",
                     "p<0.001"= "darkgreen",
                     "p<0.01" = "green",
                     "p<0.05" =  "lightgreen"), 
          title = "Nivel de significación",
          size = 0.5) +
    tm_shape(provincias) +
    tm_borders()

Y como curiosidad, hasta podríamos crear un visor con esta información

tmap_mode("view")
tmap mode set to interactive viewing
tm_shape(lisa) +
    tm_dots(col = "mean",
          style = "cat",
          palette= c("NS"="grey",
                     "Low-Low" = "blue",        # Low-Low
                     "Low-High" = "lightblue",  # Low-High 
                     "High-Low"= "sienna1",     # High-low
                     "High-High"="red"),        # High -High
          title = "LISA con p > 0.1") +
    tm_shape(provincias) +
    tm_borders(col="white") +
    tm_basemap(leaflet::providers$Esri.WorldImagery)

7.5.2 I local de Anselin con ArcGis Pro

ArcMap permite realizar el equivalente a un análisis cluster LISA de Anselin con suma facilidad. Para ello debemos seleccionar la herramienta Cluster and Outlier Analysis (Anselin Local Moran’s I) que se encuentra en el menú Spatial Statistics Tools/Mapping Clusters.

El proceso que sigue ArcMap es similar: calcula la I local, los z-scores, el p-valor, y finalmente realiza una agrupación en clusters como hemos visto antes:

Sin embargo, hay que poner atención en un aspecto: ArcMap identifica los puntos con un sistema de colores inverso al que hemos visto anteriormente. Para ArcMap, los valores High-Low y Low-High son outliers, y por lo tanto los identifica con colores más vivos, mientras que los High-High y Low-Low, que son los que se destacan normalmente en un análisis LISA, aquí se representan con tonos más pálidos.

7.5.3 I local de Anselin con GeoDa

Igual que hicimos antes, vamos a ver como realizar un análisis LISA con GeoDa, en lugar de hacerlo con QGis. Como no podía ser de otra manera este es uno de los puntos fuertes de este software, creado por el propio Luc Anselin. Para ello debemos buscar el comando Univariate Local Moran's I, en el menú Space. Seleccionando la variable de interés y el fichero de vecinos, podemos generar todos los elementos del LISA. Uno de los atractivos de GeoDa es que estos gráficos son interactivos y están conectados, de manera que seleccionar unos puntos en uno de ellos resalzará los mismos puntos en el resto de gráficos:

Por supuesto, estas no son las únicas funcionalidades de GeoDa, como veremos más adelante durante el curso.